
package es.unex.sextante.libMath.regression;

import java.text.DecimalFormat;
import java.util.ArrayList;

import Jama.Matrix;
import Jama.QRDecomposition;

/**
 * Least squares fitting
 * @author volaya
 *
 */
public class LeastSquaresFit {

    private ArrayList m_X = new ArrayList();

    private ArrayList m_Y = new ArrayList();

    private double[] m_dCoeffs;

    private String m_sExpression;

    private double m_dYMin = Double.MAX_VALUE;

    private double m_dXMax = Double.NEGATIVE_INFINITY;

    private double m_dXMin = Double.MAX_VALUE;

    private double m_dYMax = Double.NEGATIVE_INFINITY;

    public void addValue(double dX, double dY) {

        setMinMaxX(dX);
        setMinMaxY(dY);

        m_X.add(new Double(dX));
        m_Y.add(new Double(dY));

    }

    public boolean calculate(int iOrder) {

        double dX, dY;
        int nk = iOrder + 1;

        double[][] alpha = new double[nk][nk];
        double[] beta = new double[nk];
        double term = 0;

        m_dCoeffs = new double[nk];

        int iNumPoints = m_X.size();

        for (int k = 0; k < nk; k++) {

            for (int j = k; j < nk; j++) {
                term = 0.0;
                alpha[k][j] = 0.0;
                for (int i = 0; i < iNumPoints; i++) {

                    double prod1 = 1.0;
                    dX = ((Double) m_X.get(i)).doubleValue();
                    if (k > 0) {
                        for (int m = 0; m < k; m++) {
                            prod1 *= dX;
                        }
                    }
                    double prod2 = 1.0;
                    if (j > 0) {
                        for (int m = 0; m < j; m++) {
                            prod2 *= dX;
                        }
                    }

                    term = (prod1 * prod2);

                    alpha[k][j] += term;
                }
                alpha[j][k] = alpha[k][j];
            }

            for (int i = 0; i < iNumPoints; i++) {
                dX = ((Double) m_X.get(i)).doubleValue();
                dY = ((Double) m_Y.get(i)).doubleValue();
                double prod1 = 1.0;
                if (k > 0) {
                    for (int m = 0; m < k; m++) {
                        prod1 *= dX;
                    }
                }
                term = (dY * prod1);

                beta[k] += term;
            }
        }

        Matrix alpha_matrix = new Matrix(alpha);
        QRDecomposition alpha_QRD = new QRDecomposition(alpha_matrix);
        Matrix beta_matrix = new Matrix(beta, nk);
        Matrix param_matrix;
        try {
            param_matrix = alpha_QRD.solve(beta_matrix);
        } catch (Exception e) {
            return false;
        }

        DecimalFormat df = new DecimalFormat("####.#####");
        StringBuffer sb = new StringBuffer("");
        for (int k = 0; k < nk; k++) {
            m_dCoeffs[k] = param_matrix.get(k, 0);
            if (k != 0) {
                sb.append(" + " + df.format(m_dCoeffs[k]) + "x^"
                        + Integer.toString(k));
            } else {
                sb.append(df.format(m_dCoeffs[k]));
            }
        }

        m_sExpression = sb.toString();

        return true;

    }

    public String getExpression() {

        return m_sExpression;

    }

    public int getNumPoints() {

        return m_X.size();
    }

    public void getPoints(double[] x, double[] y) {

        int i;

        for (i = 0; i < m_X.size(); i++) {
            x[i] = ((Double) m_X.get(i)).doubleValue();
            y[i] = ((Double) m_Y.get(i)).doubleValue();
        }

    }

    public double getXMax() {

        return m_dXMax;

    }

    public double getXMin() {

        return m_dXMin;

    }

    public double getY(double x) {

        int i;
        double dRet = 0;

        for (i = 0; i < m_dCoeffs.length; i++) {
            dRet += m_dCoeffs[i] * Math.pow(x, (double) i);
        }

        return dRet;

    }

    public double getYMax() {

        return m_dYMax;

    }

    public double getYMin() {

        return m_dYMin;

    }

    private void setMinMaxX(double x) {

        if (x > m_dXMax) {
            m_dXMax = x;
        }
        if (x < m_dXMin) {
            m_dXMin = x;
        }

    }

    private void setMinMaxY(double y) {

        if (y > m_dYMax) {
            m_dYMax = y;
        }
        if (y < m_dYMin) {
            m_dYMin = y;
        }

    }
}
